{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "10bf9962-b18d-47a5-80d4-b66348dd8956",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# The paleopiezometry module\n",
    "\n",
    "**What is paleopizometry?**  \n",
    "Paleopiezometers refer to microstructural features of deformed rocks that vary with the magnitude of the applied differential stress under which they formed (Twiss and Moores, 2006). They serve to infer differential stresses in the geological past, hence the name paleopiezometry, and they are an essential tool for validating rheological models of the lithosphere. The most prevalent microstructure used for this purpose is the average recrystallized (apparent) grain size. This choice is due to its ease of measurement in recrystallized rocks.\n",
    "\n",
    "> 📣 The GrainSizeTools script includes a **Jupyter notebook template** to promote the reproducibility of palaeopizometry studies. Use it as a template by deleting and adding as necessary and use it as supplementary material to your study so that anyone can reproduce your results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "93287c25-1187-44c0-985d-3f293cf99762",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# modules import\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e05b390",
   "metadata": {},
   "source": [
    "## Estimate the stress using the average grain size\n",
    "\n",
    "### Before you begin\n",
    "\n",
    "The paleopiezometry module includes a function named ``calc_diffstress()`` for estimating differential stress based on average recrystallized grain sizes. The function requires careful consideration of the following factors for correct differential stress estimation:\n",
    "\n",
    "1. definition of the mineral phase and the piezometer relationship to be used\n",
    "2. enter the (apparent) grain size as **equivalent circular diameters** in microns\n",
    "3. enter a specific average (arithmetic mean, median, etc.) with **no stereological correction**,\n",
    "4. set the stress framework, either uniaxial compression/extension or plane stress, for proper stress correction.\n",
    "\n",
    "Regarding the first factor, GrainSizeTools includes a database of paleopiezometric relations, including the most common mineral phases quartz, calcite, olivine and albite feldspar. In addition, the script facilitates the writing of ad hoc piezometric relations. As we will see later, it is also very easy to add new piezometric relationships.\n",
    "\n",
    "For the second factor, the function assumes by default that the average grain size you specify has been estimated using equivalent circular diameters (ECDs). As some piezometric relationships, especially those established longer ago, were established using linear intercepts (LI) rather than ECDs, the function will automatically convert ECDs to linear intercepts if necessary using the Hoff and Rhines (1968) correction. This means that **you do not need to worry about whether the piezometer was originally calibrated using linear intercepts or not**, always use averages based on equivalent circular diameters in microns. The function will explicitly warn you when using this ECD ti LI of conversion.\n",
    "\n",
    "The third factor is key for the correct estimation of differential stress, as each paleopiezometry relationship has been calibrated to a specific average grain size (e.g. the arithmetic mean, median or RMS mean) and thus **only gives valid results if the same type of average is used**. In addition, **you must not apply any type of stereological correction to the calculated grain average**, if the author(s) of the piezometer used any type of stereological correction during calibration, your input will be automatically corrected by the function.\n",
    "\n",
    "The fourth factor means that the user must decide whether or not to correct the differential stress estimate for plane strain using the correction factor suggested by Paterson and Olgaard (2000). The reason for this is that piezometer calibration experiments are mainly performed in uniaxial compression, whereas natural shear zones behave approximately like plane-strain volumes.\n",
    "\n",
    "### Use of the piezometric database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7fa91a85",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "module plot imported\n",
      "module averages imported\n",
      "module stereology imported\n",
      "module piezometers imported\n",
      "module template imported\n",
      "\n",
      "======================================================================================\n",
      "Welcome to GrainSizeTools script\n",
      "======================================================================================\n",
      "A free open-source cross-platform script to visualize and characterize grain size\n",
      "population and estimate differential stress via paleopizometers.\n",
      "\n",
      "Version: 2024.02.xx\n",
      "Documentation: https://marcoalopez.github.io/GrainSizeTools/\n",
      "\n",
      "Type get.functions_list() to get a list of the main methods\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Load the script first (change the path to GrainSizeTools_script.py accordingly!)\n",
    "%run C:/Users/marco/Documents/GitHub/GrainSizeTools/grain_size_tools/GrainSizeTools_script.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "630098b8",
   "metadata": {},
   "source": [
    "You can get information from the console on the different available piezometric relations  just by typing ``piezometers._()``, where \\_ is the mineral phase, either ``quartz``, ``calcite``, ``olivine``, or ``feldspar``. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f5787f0b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Available piezometers:\n",
      "'Cross'\n",
      "'Cross_hr'\n",
      "'Holyoke'\n",
      "'Holyoke_BLG'\n",
      "'Shimizu'\n",
      "'Stipp_Tullis'\n",
      "'Stipp_Tullis_BLG'\n",
      "'Twiss'\n"
     ]
    }
   ],
   "source": [
    "piezometers.quartz()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f33f79e6",
   "metadata": {},
   "source": [
    "If you want to get the details of a particular piezometric relationship simply pass the name of the relationship inside the parenthesis. Remember that the relationship between recrystallized grain size and differential stress is\n",
    "\n",
    "$$\\sigma_d = B\\ d^{-m}$$\n",
    "\n",
    "where $\\sigma_d$ and $d$ are the differential stress and the average grain size respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ec3b87b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(550,\n",
       " 0.68,\n",
       " 'Ensure that you entered the apparent grain size as the arithmetic mean grain size',\n",
       " True,\n",
       " 1.5)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "piezometers.quartz('Twiss')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d478c6ff",
   "metadata": {},
   "source": [
    "The five different outputs separated by commas correspond with:\n",
    "\n",
    "- the constant *B* of the piezometric relation\n",
    "- the exponent *m* of the piezometric relation\n",
    "- A warning indicating the average to use with this piezometric relation\n",
    "- An indication of whether the piezometric relation was calibrated using linear intercepts (if ``False`` the piezometric relation was calibrated using equivalent circular diameters).\n",
    "- The stereological correction factor used (if applicable). If ``False``, no stereological correction applies.\n",
    "\n",
    "### Estimating differential stresses\n",
    "\n",
    "Let us first look at the documentation of the:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "672836d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1;31mSignature:\u001b[0m \u001b[0mcalc_diffstress\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgrain_size\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mphase\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpiezometer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcorrection\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mDocstring:\u001b[0m\n",
      "Apply different piezometric relations to estimate the differential\n",
      "stress from average apparent grain sizes. The piezometric relation has\n",
      "the following general form:\n",
      "\n",
      "df = B * grain_size**-m\n",
      "\n",
      "where df is the differential stress in [MPa], B is an experimentally\n",
      "derived parameter in [MPa micron**m], grain_size is the aparent grain\n",
      "size in [microns], and m is an experimentally derived exponent.\n",
      "\n",
      "Parameters\n",
      "----------\n",
      "grain_size : positive scalar or array-like\n",
      "    the apparent grain size in microns\n",
      "\n",
      "phase : string {'quartz', 'olivine', 'calcite', or 'feldspar'}\n",
      "    the mineral phase\n",
      "\n",
      "piezometer : string\n",
      "    the piezometric relation\n",
      "\n",
      "correction : bool, default False\n",
      "    correct the stress values for plane stress (Paterson and Olgaard, 2000)\n",
      "\n",
      " References\n",
      "-----------\n",
      "Paterson and Olgaard (2000) https://doi.org/10.1016/S0191-8141(00)00042-0\n",
      "de Hoff and Rhines (1968) Quantitative Microscopy. Mcgraw-Hill. New York.\n",
      "\n",
      "Call functions\n",
      "--------------\n",
      "piezometers.quartz\n",
      "piezometers.olivine\n",
      "piezometers.calcite\n",
      "piezometers.albite\n",
      "\n",
      "Assumptions\n",
      "-----------\n",
      "- Independence of temperature (excepting Shimizu piezometer), total strain,\n",
      "flow stress, and water content.\n",
      "- Recrystallized grains are equidimensional or close to equidimensional when\n",
      "using a single section.\n",
      "- The piezometer relations requires entering the grain size as \"average\"\n",
      "apparent grain size values calculated using equivalent circular diameters\n",
      "(ECD) with no stereological correction. See documentation for more details.\n",
      "- When required, the grain size value will be converted from ECD to linear\n",
      "intercept (LI) using a correction factor based on de Hoff and Rhines (1968):\n",
      "LI = (correction factor / sqrt(4/pi)) * ECD\n",
      "- Stress estimates can be corrected from uniaxial compression (experiments)\n",
      "to plane strain (nature) multiplying the paleopiezometer by 2/sqrt(3)\n",
      "(Paterson and Olgaard, 2000)\n",
      "\n",
      "Returns\n",
      "-------\n",
      "The differential stress in MPa (a float)\n",
      "\u001b[1;31mFile:\u001b[0m      c:\\users\\marco\\documents\\github\\grainsizetools\\grain_size_tools\\grainsizetools_script.py\n",
      "\u001b[1;31mType:\u001b[0m      function"
     ]
    }
   ],
   "source": [
    "calc_diffstress?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a2db1fe",
   "metadata": {},
   "source": [
    "As indicated in the documentation, the ``calc_diffstress()`` requires at least three inputs: (1) the average grain size in microns, (2) the mineral phase, and (3) the piezometric relation to be used. A few examples are given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e5f88969",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "============================================================================\n",
      "differential stress = 83.65 MPa\n",
      "\n",
      "INFO:\n",
      "Ensure that you entered the apparent grain size as the arithmetic mean grain size\n",
      "ECD was converted to linear intercepts using de Hoff and Rhines (1968) correction\n",
      "============================================================================\n"
     ]
    }
   ],
   "source": [
    "calc_diffstress(12, phase='quartz', piezometer='Twiss')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6806817",
   "metadata": {},
   "source": [
    "The function returns the calculated differential stress (in MPa) plus some relevant information about the corrections made and the type of average expected as input.\n",
    "\n",
    "As most piezometric calibrations have been calibrated using uniaxial compression deformation experiments, let's correct this estimate for planar stress using the correction suggested by Paterson and Olgaard (2000). This is done by passing a new parameter within the function as follows (note the slightly different value of differential stress):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ee74c503",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "============================================================================\n",
      "differential stress = 96.59 MPa\n",
      "\n",
      "INFO:\n",
      "Ensure that you entered the apparent grain size as the arithmetic mean grain size\n",
      "ECD was converted to linear intercepts using de Hoff and Rhines (1968) correction\n",
      "============================================================================\n"
     ]
    }
   ],
   "source": [
    "calc_diffstress(12, phase='quartz', piezometer='Twiss', correction=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efc15b70",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## How to synthesise a set of piezometric values\n",
    "\n",
    "TODO\n",
    "\n",
    "> ⚠️ **Never combine averages -or any other statistic- as if they were data points**. In the context of meta-analysis, which is a statistical procedure for combining data from multiple studies or experiments, a common mistake is to take different estimates and averaged them. However, the stress estimates are already an average and calculating a mean of the averages is mathematically incorrect because this approach ignores the uncertainty of the average. Similarly, to look at the raw grain size data from all the grain maps and make the estimate from there is also incorrect because this ignores the uncertainty within and between experimental conditions. The correct procedure here is to take the weighted mean of the averages, where each mean is weighted by its variance or the squared standard error of the mean (SEM). The GrainSize Tools script has a function to do this named ```weighted_mean_and_se()```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a697822",
   "metadata": {},
   "source": [
    "## Summary of piezometric relationships in the database\n",
    "\n",
    "$$\\sigma_d = B\\ d^{-m}$$\n",
    "$$d = A\\ \\sigma_d^{-p}$$\n",
    "\n",
    "$B$ and $m$ relate to $A$ and $p$ as follows:\n",
    " \n",
    "$$B = A^{1/p}$$\n",
    "$$m = 1/p$$\n",
    "\n",
    "$A$ and $B$ are in $\\mu m\\ MPa^{p,m}$\n",
    "\n",
    "### Quartz piezometric relationships\n",
    "\n",
    "|          calibration           |       reference        |    DRX    |        average         |    A    |  p   |   B    |  m   |\n",
    "| :----------------------------: | :--------------------: | :-------: | :--------------------: | :-----: | :--: | :----: | :--: |\n",
    "|     Cross et al. (2017)*§*     |      ``'Cross'``       | BLG, SGR  |        RMS mean        | 8128.3  | 1.41 | 593.0  | 0.71 |\n",
    "|     Cross et al. (2017)*§*     |     ``'Cross_hr'``     | BLG, SGR  |        RMS mean        | 16595.9 | 1.59 | 450.9  | 0.63 |\n",
    "| Holyoke & Kronenberg (2010)*¶* |     ``'Holyoke'``      | SGR + GBM |        RMS mean        |  2451   | 1.26 | 490.3  | 0.79 |\n",
    "| Holyoke & Kronenberg (2010)*¶* |   ``'Holyoke_BLG'``    |    BLG    |        RMS mean        |   39    | 0.54 | 883.9  | 1.85 |\n",
    "|        Shimizu (2008)‡         |     ``'Shimizu'``      | SGR + GBM | Median in log(e) scale |  1525   | 1.25 |  352   | 0.8  |\n",
    "|    Stipp and Tullis (2003)     |   ``'Stipp_Tullis'``   | SGR + GBM |        RMS mean        | 3630.8  | 1.26 | 669.0  | 0.79 |\n",
    "|    Stipp and Tullis (2003)     | ``'Stipp_Tullis_BLG'`` |    BLG    |        RMS mean        |   78    | 0.61 | 1264.1 | 1.64 |\n",
    "|         Twiss (1977)†          |      ``'Twiss'``       | SGR + GBM |      arith. mean       |  1230   | 1.47 |  550   | 0.68 |\n",
    "\n",
    "*† Apparent grain size measured as equivalent circular diameters (ECD) with no stereological correction and reported in microns. The use of non-linear scales is indicated*    \n",
    "*‡ Shimizu piezometer requires to provide the temperature during deformation in K*    \n",
    "*§ Cross et al. (2017) reanalysed the samples of Stipp and Tullis (2003) using EBSD data for reconstructing the grains. Specifically, they use grain maps with a 1 m and a 200 nm (hr - high-resolution) step sizes . This is the preferred piezometer for quartz when grain size data comes from EBSD maps*  \n",
    "*¶ Holyoke and Kronenberg (2010) provides a recalibration of the Stipp and Tullis (2003) piezometer*  \n",
    "\n",
    "> To check:  \n",
    "> Tockle and Hirth 2021 https://doi.org/10.1029/2020JB021475  \n",
    "> Bishop (1996)  \n",
    "> Kidder et al. (2016) https://doi.org/10.1016/j.jsg.2015.12.004  \n",
    "> Heilbronner & Kilian (2017)  \n",
    "> Richter et al. (2018)  \n",
    "> Soleymani et al. (2020) https://doi.org/10.1130/G46972.1\n",
    "\n",
    "### Olivine piezometric relationships\n",
    "\n",
    "|        calibration         |      reference      | DRX  | dry/wet |   average   |  A    |  p   |    B    |  m   |\n",
    "| :------------------------: | :-----------------: | :--: | :-----: | :---------: | :---: | :--: | :-----: | :--: |\n",
    "|  Jung and Karato (2001)§   |  ``'Jung_Karato'``  | BLG  |   wet   | arith. mean | 25704 | 1.18 | 5461.03 | 0.85 |\n",
    "| Van der Wal et al. (1993)§ | ``'VanderWal_wet'`` |      | dry/wet | arith. mean | 15000 | 1.33 | 1355.4  | 0.75 |\n",
    "|    Tasaka et al. (2015)    |  ``'Tasaka_wet'``   |      |   wet   | arith. mean | 6310  | 1.33 |  719.7  | 0.75 |\n",
    "\n",
    "*§ These piezometers were originally calibrated using linear intercepts (LI) instead of ECD*  \n",
    "\n",
    "> To check: Add Twiss 1977, Karato 1980 and Ross et al 1980?\n",
    "\n",
    "\n",
    "### Calcite piezometric relationships\n",
    "\n",
    "|        calibration          |      reference      |   DRX    |   average   |   A    |  p   |    B    |  m   |\n",
    "| :-------------------------: | :-----------------: | :------: | :---------: | :----: | :--: | :-----: | :--: |\n",
    "|   Barnhoorn et al. (2004)   |   ``'Barnhoorn'``   | SRG, GBM | arith. mean | 2134.4 | 1.22 | 537.03  | 0.82 |\n",
    "| Platt and De Bresser (2017) | ``'Platt_Bresser'`` | BLG, SGR |  RMS mean   |  2141  | 1.22 | 538.40  | 0.82 |\n",
    "|        Rutter (1995)        |  ``'Rutter_SGR'``   |   SGR    | arith. mean | 2026.8 | 1.14 | 812.83  | 0.88 |\n",
    "|        Rutter (1995)        |  ``'Rutter_GBM'``   |   GBM    | arith. mean | 7143.8 | 1.12 | 2691.53 | 0.89 |\n",
    "|    Valcke et al. (2015)     |    ``'Valcke'``     | BLG, SGR | arith. mean | 79.43  | 0.6  | 1467.92 | 1.67 |  \n",
    "\n",
    "> To check:: Add Twiss 1977 and Schmidt et al. 1980?\n",
    "\n",
    "### Feldspar piezometric relationships\n",
    "\n",
    "|       calibration       |       reference       | DRX  | average |  A   |  p   |   B   |  m   |\n",
    "| :---------------------: | :-------------------: | :--: | :-----: | :--: | :--: | :---: | :--: |\n",
    "| Post and Tullis (1999)§ | ``'Post_Tullis_BLG'`` | BLG  | Median  |  55  | 0.66 | 433.4 | 1.52 |\n",
    "\n",
    "*§ These piezometers were originally calibrated using linear intercepts (LI) instead of ECD*\n",
    "\n",
    "> TODO:  \n",
    "> Add Speciale et al. (2022) https://doi.org/10.1016/j.jsg.2021.104495\n",
    "\n",
    "\n",
    "### Orthopyroxene piezometric relationships\n",
    "\n",
    "> TODO:  \n",
    "> Linckens et al. 2014 https://doi.org/10.1016/j.epsl.2013.11.037  \n",
    "> Brujin ans Skemer 2014 https://doi.org/10.1002/2014GL060607  \n",
    "> Speciale et al. 2022 https://doi.org/10.1016/j.jsg.2021.104495\n",
    "\n",
    "\n",
    "### Halite piezometric relationships\n",
    "\n",
    "> TODO:  \n",
    "> Guillopé and Poirier 1979, https://doi.org/10.1029/JB084iB10p05557  \n",
    "> Ter Heege et al. 2005 https://doi.org/10.1016/j.tecto.2004.10.002\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d54696e5-6445-4351-833d-8f1b9a865313",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook tested in 2023-12-22 using:\n",
      "Python 3.10.13 | packaged by Anaconda, Inc. | (main, Sep 11 2023, 13:15:57) [MSC v.1916 64 bit (AMD64)]\n",
      "Numpy 1.26.2\n",
      "Matplotlib 3.8.0\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "from datetime import date    \n",
    "today = date.today().isoformat()\n",
    "import matplotlib as mpl\n",
    "\n",
    "print(f'Notebook tested in {today} using:')\n",
    "print('Python', sys.version)\n",
    "print('Numpy', np.__version__)\n",
    "print('Matplotlib', mpl.__version__)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
